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1.   Introduction 

In  a  recent  report,  Zlatev  [4]  extends  the  known  class  of 

unconditionally  stable  variable-step  multistep  methods  (that  is,  methods 

which  are  stable  for  any  stepsize  sequence  such  that  the  ratio  of  the 

largest  to  smallest  steps  is  bounded)  to  include  methods  of  the  form 

(1.1)      y  =  0L  [{h}]y  .  4-  a.[{h}]y   „  +  derivative  terms, 
n    i      n-1    2.  n-Z 

By  the  notation  a.[{h}]  we  simply  indicate  that  the  a.  may  depend  on  the 

ratio  of  one  or  more  past  steps.   Zlatev  shows  (1.1)  to  be  stable  as  long  as 

|  Oi,,  |  <  1,  a  property  that  will  hold  for  any  reasonable  strategy  for 

computing  a,  and  a_  from  {h}.   (Note  that  |a„|  <  1  for  a  strongly  stable 

constant-step  method.) 

The  proof  of  convergence  presented  in  Zlatev' s  report  is 
particularly  interesting  for  two  reasons:   (a)  it  yields  a  proof  of 
fixed-step  convergence  that  is  intuitively  more  pleasing  than  the 
"standard"  proof,  as  found  in  Henrici  [3],  being  more  direct;  and 
(b)  it  permits  an  extension  to  all  strongly  stable  methods,  which,  while 
not  yielding  unconditional  stability  (shown  to  be  untrue  in  general  by  a 
counter  example  in  [1]),  leads  to  a  bound  on  stepsize  ratios  which  can 
be  computed  for  any  particular  formula,  and  hence  simply  implemented  in 
practical  programs. 

This  short  paper  discusses  an  extension  of  Zlatev' s  proof, 
first  at  an  intuitive  level  which  this  writer  finds  particularly 
appealing,  then  at  a  technical  level.   Finally,  the  constructive 
application  of  the  result  to  the  calculation  of  stable-step  ratios 
(the  bounds  b  and  B)  is  discussed. 


2.   Zlatev's  Proof  Technique 

Consider  the  strongly  stable,  constant-step  method 


(2.1)  y  =  Z  a.  y   .  +  Z  h3.  f(y   .) 

n    .  ..   1  n-i  .  _   i    n-i 

i=l  i=0 

for  the  initial  value  problem 

(2.2)  y'  =  f(y,  t),  y(0)  =  yQ. 

The  equation  for  the  global  error,  e  =  y  -  y(t  ),  takes  the  form 

n    n      n 


(2.3)  e  =  Z  a.e   .+  I     h3.  g   .e   .  +  d   , 

n    .  ,   l  n-i   .  _   i  n-i  n-i    n 
i=l  i=0 

where  g.  =  f  evaluated  at  a  suitable  point,  and  d   is  the  local  truncation 
:    y  n 

error.   The  "standard"  proof  of  convergence  involves  bounding  the  rate  of 

growth  of  solutions  of  this  difference  equation  which  is  an  order  h 

perturbation  of  a  difference  equation  with  a  root  at  one.   Zlatev  considers, 

instead,  a  difference  equation  for  c  =  e  -  e   „  .   This  satisfies 

n    n    n-1 

the  equation 

k-1  k 

(2.4)  c  =e  -e  ,=  Z  y.    c   .+  Z  hS .  g   .e   ,  +  d 

n    n    n-1    .    '  i  n-i    .  „   l  n-i  n-i    n 
i=l  i=0 

If  we  write  p(£)  =  -£  +  £  & .    K        ,   P(£)  =  ~K     +     Z  Y.C    we  find 

i=l  1  i=l   1 

that  p(£)  =  K   P(?)/(?-D,  that  is,  p(?)  is  a  "deflated"  form  of  p(0 

with  the  unit  root  replaced  by  zero.   If  the  method  is  strongly  stable, 

p(£)  has  all  of  its  zeros  strictly  inside  the  unit  circle,  hence  the 

growth  of  solutions  of  eq.  (2.4)  can  be  bounded  by  a  constant  times  the 

__  i -i 

largest  of  the  added  terms  (  Z  h$ .  g   .  e   .  +  d  ).   If  d  is  0(h   ) 

.  _   i  n-i  n-i    n        n 
i=0 

and  e  is  0(hP),  this  term  is  0(hP+1),  hence  c   is  0(hP   ).   The  global 


n 

error  e  is  simply  en  +     Z  (c.),  and,  since  there  are  order  of  1/h  terms 
n  U    .  t   1 

i=l 

in  the  sum,  e  is  0(h   /h)  =  0(h).   Of  course,  this  argument  is 

circular,  so  there  are  a  few  technical  details  to  work  out.   This  is 

the  subject  of  the  next  section. 


3.   Variable-step  Methods 

Consider  a  method  given  by  eq.  (2.1)  where  h  depends  on  n,  and 
the  a.  and  $.  may  vary  with  the  ratios  of  stepsizes.   Following  the 
notation  in  [1],  the  error  equation  (2.3)  may  be  written  in  the  form 
(3.1) 


e  =(S  +hS)e   n+d 
— n     n    n  n  — n-1   — n 


where  e  is  the  error  in  the  collection  of  saved  information 


*n  (e'g' 


^n  =  [V  yn- 


1'  V-2 


, .  . .  ,  h  y '  ,  h 


nyn'   n-1  'n-l 


, .  .  .  ]  )  ,  S   contains  the 


dependence  of  (2.3)  on  g  =  3f/8y,  and  S  depends  only  on  the  coefficients 


of  the  method.   In  [1]  it  is  shown  that  if   S    (where 

m 

S  , n )  and   I S    can  be  bounded  independently 

m+1       ' '  n' ' 


n 


sn  =     s  =s  s  ,. 

m   .   , _   1   n  n-1 

j=m+l 


of  n  >  m,  then  convergence  can  be  proved  for  consistent  methods,  so  that 
bounded   | S    and   | S    is  a  reasonable  stability  condition.   For  the 


m       '  '  n 
multistep  method  above, 


s  = 
n 


\    a2     ...           akj31     62     ... 

K 

10...              0 

0       1...              0                 ( ) 

0        0                  10 

_                      _      1 

r 

0        0 

0 

1        0       ... 

O    |°  — 

0 
0 

• 

0        0 

1       0_ 

We  now  consider  the  similarity  transformation  of  S   to  Q  S  Q   where 

n       n 


Q  - 


l  o  ... 
-l  l  ... 


O 


-l    i 


o 


o 


This  corresponds  to  a  change  of  variables,  from  e,e  ,,e  O,...to 

n   n-1   n-z 

e  ,  -c  ,  -c  ,,...  and  is  in  the  spirit  of  Zlatev's  technique.   We  find  that 
n    n    n-1 


QSnQ 


-1 


1 

"Yl     "Y2    ••'   "Yk-1 

3,      3«     ... 

\ 

0 

Yl       Y2    •'•     Yk-1 

-3 ,     -3  _     ... 

A 

0 

10                   0 
1          0 

L 

O 

1 

0          0        ... 

0 

o 

1          0 
0          1 

0 
0 

1        0_ 

n 


n 


sn  q  1  =   n   (Q  s.  q"1) 

111  Li  J 

j=m+l 


1  !        y              ! 

n,m 
— t , 

^n                     1 
0   I                 S                                          B 

m                                        n,m 

C^)      \      cn~m 

where 
(3.2) 


y    =  y  .    +  x  S     with  y 
n ,  m    n—  1 ,  m    n  m         mm 


=  0  . 


Suppose,  for  a  moment,  that   S    <  r  where  r  <  1. 

n   — 


Then  eq.  (3.2)  implies 


(3.3) 


n,m' '  —  1-r 


max 
m<i<n 


Thus,  if  the  y .    are  bounded,  the  matrix 


i  y 
o 


n,m 

n 

m 


,n-m 


is  bounded.   Because  C    =0  for  n  >  m  +  k  and  is  trivially  bounded,  it 

then  follows  that  B    is  bounded  if  the  3.  are  bounded.   Thus,  it  remains 

n,m  x 

to  consider  conditions  under  which   Is    <  r  <  1.   Note  that  S   is  the 

ii  ni  i  _  n 

companion  matrix  of  the  polynomial  p(£)/(£-l),  all  °f  whose  zeros  are 
inside  the  unit  circle.  Consider  the  coefficients  for  eq.  (2.1)  with 
constant  steps.   If  the  largest  zero  of  p (£)/(?-!)  is  r,  there  exists  a 


norm  | | • | |   such  that 


S| |   <  r  +  E  =  r 


for  any  £  >  0  where  S  is  S   for  constant  stepsizes.   (Let  H  be  a  similarity 

transformation  S  to  Jordan  form  with  its  nonzero  off-diagonal  elements 

equal  to  e.   Let   |x| I   be  the  max  norm  of   Hx   and   si  I   the  matrix  norm 

n  H 

I.  /\    .  . 
S   TT.   If  the 
1  n '  '  H 

coefficients  a.  of  a   variable-step  method  are  continuous  functions  of  the 
1 

|S     is  a  continuous 
n  H 

function  of  the  step  ratios.   When  all  the  ratios  are  unity, 

|§  | |   <  r  <  1.   Hence,  for  an  r  such  that  r  <  r  <  1,  there  exists  a  pair 
n  H 

of  numbers  b  <  1  <  B  such  that  if  the  step  ratios  lie  between  b  and  B 

|S    <  r  <  1. 
ii  u i i  _ 


4.   Computing  b  and  B 

This  can  be  done  once  and  for  all  for  any  given  method,  it 
does  not  have  to  be  computed  during  an  integration.   First,  the  similarity 

transform  H  of  S  to  Jordan  form  must  be  computed.   (If  the  zeros  of  p(£) 

i  i  ^  ii 
are  distinct,  the  Jordan  form  is  diagonal.)   Next,  the  function   |S  | |„, 

which  depends  on  the  step  ratios,  must  be  formed.   For  any  pair  of  values 

b  and  B.  the  maximum  of   S  I  LT  in  the  hypercube  of  step  ratios  between  b 
/  ' '  n' 'H 

and  B  can  be  computed  (approximately) .   The  value  of  b  and  B  can 

be  decreased  and  increased  respectively  until  this  maximum 

is  any  selected  r  such  that  1  >  r  >  f.   (If  r  is  chosen  closer  to  1,  the 

range  between  b  and  B  will  be  larger,  but  the  error  constants  will  be 

greater  because  of  the  1/(1  -  r)  factor  in  eq.  (3.3).)   For  many  methods 

one  will  find  that  b  can  be  made  very  close  to  zero  without  problem  as 

rapid  step  reductions  in  variable  step  methods  yield  methods  that  are 


essentially  one-step  methods  of  order  one  or  two  for  which  the  y.  are 
zero.  This  is  fortunate  because  it  is  desirable  to  be  able  to  reduce 
the  step  suddenly  when  the  solution  exhibits  a  sharp  change. 

Variable  Order  Methods 

In  [2]  it  was  shown  that  changing  formulas  "occassionally" 

did  not  affect  stability.   Precisely,  it  was  shown  that,  for  any  set  of 

strongly  stable  formulas  {F. }  there  exists  integers  N. .  such  that  if  at 

1  ji 

least  N.  .  steps  of  formula  F.  are  taken  after  a  switch  from  formula  F., 
ji  i  J 

then  the  combined  method  is  stable.   This  result  carries  over  to  the 

Ai 
current  situation  with  a  much  simpler  proof.   Suppose  S   is  the  matrix 

S  corresponding  to  formula  F.  at  constant  stepsize.   Let  H.  be  such  that 

Is1!!    <  r.  <  1.   Choose  the  r.,  b.  and  B.  such  that   Is1!!    <  r.  <  1 
1  '   '  'H.  —  i  11      i  M  n"H,  -  i 

i  i 

for  all  step  ratios  in  the  range  [b.,  B.].   Now  choose  r  such  that 

max  r.  <  r  <  1.   (We  only  consider  finite  sets  of  formula  in  practice. 

i 

The  mathematically  minded  must  add  requirements  such  as  uniformly 

strongly  stable  to  guarantee  that  sup  r.  <  1.)   Now  consider  the  max  norm 

i 

|#|   of  S  where  all  steps  are  taken  with  method  i. 
■  i  i  i     m 

I  |cn|  I  <  I  liT1!  I   I  |qn!  I     I  |W  I  I 

i i  mi  i  _  i i  x    i  i  i  '  m'  'h.  ' '  i'  ' 

l 

MhT1!!  Mh.II  rn_m 

—  '  '  i  '  '  '  '  1' '   i 

<  UIh:1!!  ||h.||  A)n"m]  rn_m 

—  '  ■  i  '  '  '  '  i1 '   r 

Since  we  have  chosen  r.  <  r,  there  exists  an  integer  N.  such  that 

i  i 

_i  r  • 

l|H.  I   I  Ih.I   (— )P  <  1    for  all  p  >  N 


10 


Thus,  if  at  least  N.  steps  are  taken  with  formula  F.,  we  can  bound 

1  1 

Is    for  a  sequence  of  variable-step,  variable-order  formulas  by 

m 

Kr    where  K  =  max(||H.  I    Ih.II).   This  can  easily  seem  to  be 

1  '  i  '  '  '  '  l1  ' 

i  i ^n i  i 
sufficient  to  show  that   S    is  bounded  as  before.   (The  proof  above 

1  '  m1  ' 

yields  N..  =  N.  for  all  j .   In  practice,  a  smaller  value  of  N..  can  be 
7      ji    i  Ji 


obtained  by  setting 


N.  .  = 
ji 


log  llHj  h,1! 


log  r  -  log  rt 


and 


K-  maxCllH"1!!  ||H  ||). 
i,j 

Only  the  presentation  of  the  proof  is  complicated,  not  the  principle.) 
Consequently,  a  practical  step  and  order  control  restriction 
which  is  guaranteed  to  be  stable  for  strongly  stable  methods  is  as  follows: 

1.  If  the  order  was  changed  within  too  few  steps,  stay  at 
the  current  order;  otherwise,  consider  changing  the  order 
in  the  next  step. 

2.  Compute  the  "optimal"  step  and  order  for  the  next  step 
using  your  favorite  algorithm. 

3.  If  the  step  increase  is  greater  than  permitted,  use 
the  maximum  allowed. 

4.  If  the  step  decrease  is  larger  than  allowed,  change  to 
a  one-step  method  to  "restart";  otherwise,  take  the 
recommended  step. 

(The  question  of  restarting  is  the  subject  of  a  paper  in  preparation  by 
this  writer.  It  is  possible  to  restart  efficiently  at  high  order  using 
a  one-step  method  which  is  automatically  stable.) 


11 


BIBLIOGRAPHY 


1.  Gear,  C.  W.  and  K.  W.  Tu,  "The  effect  of  variable  mesh  size  on  the 

stability  of  multistep  methods,"  SIAM  Journal  on  Numerical 
Analysis,  11  (5),  October  1974,  1025-1043. 

2.  Gear,  C.  W.  and  D.  S.  Watanabe,  "Stability  and  convergence  of 

variable  order  multistep  methods,"  SIAM  Journal  on  Numerical 
Analysis,  11  (5),  October  1974,  1044-1058. 

3.  Henrici,  P.,  DISCRETE  VARIABLE  METHODS  IN  ORDINARY  DIFFERENTIAL 

EQUATIONS,  John  Wiley  &  Sons,  Inc.,  New  York,  1962. 

4.  Zlatev,  Z.,  "Stability  properties  of  variable  step-size,  variable 

formula  methods,"  Mathematics  Institute,  University  of 
Copenhagen,  Preprint  Series  //38,  1977. 


„mAEC-427  U.S.  ATOMIC   ENERGY  COMMISSION 

SfSm  UNIVERSITY-TYPE  CONTRACTOR'S  RECOMMENDATION   FOR 

DISPOSITION  OF  SCIENTIF  C  AND  TECHNICAL  DOCUMENT 

I  See  Instructions  on  Reverse  Side  ) 


AEC  REPORT  NO. 

COO-2383-0048 


STABILITY  OF  VARIABLE-STEP  METHODS 
FOR  ORDINARY  DIFFERENTIAL  EQUATIONS 


TYPE  OF   DOCUMENT    (Check  one): 

[x]  a.  Scientific  and  technical  report 
[]  b.  Conference  paper  not  to  be  published  in  a  journal: 

Title  of  conference 

Date  of  conference 


Exact  location  of  conference _ 

Sponsoring  organization 

□  c.  Other   (Specify) 


RECOMMENDED  ANNOUNCEMENT  AND  DISTRIBUTION   (Check  one): 

Kl  a.  AEC's  normal  announcement  and  distribution  procedures  may  be  followed. 

J  b.   Make  available  only  within  AEC  and  to  AEC  contractors  and  other  U.S.  Government  agencies  and  their  contractors. 
]  c.   Make  no  announcement  or  distribution. 

REASON    FOR    RECOMMENDED    RESTRICTIONS: 


SUBMITTED   BY:      NAME   AND  POSITION   (Please  print  or  type) 

C.   William  Gear 
Principal   Investigator 


Organization 

Department   of   Computer   Science 

University  of   Illinois   at  Urbana-Champaign 

Urbana,    Illinois   61801 


Signature 


Date 

July  18,    1978 


FOR   AEC   USE   ONLY 

AEC  CONTRACT  ADMINISTRATOR'S  COMMENTS,   IF    ANY,  ON    ABOVE   ANNOUNCEMENT  AND   DISTRIBUTION 
RECOMMENDATION: 


PATENT  CLEARANCE: 

Lj  a.  AEC  patent  clearance  has  been  granted  by  responsible  AEC  patent  group. 
LJ  b.   Report  has  been  sent  to  responsible  AEC  patent  group  for  clearance. 
LJ  c.  Patent  clearance  not  required. 


BLIOGRAPHIC  DATA 
IEET 

1.   Report  No. 

UIUCDCS-R-78-931 

2. 

3.  Recipient's  Accession  No. 

Title  and  Subtitle 

STABILITY   OF   VARIABLE-STEP  METHODS 

5.  Report  Date 

Julv  18,    1978 

FOR  ORDINARY   DIFFERENTIAL   EQUATIONS 

6. 

Author(s) 

C.    William  Gear 

8.   Performing  Organization  Rept. 

No-UIUCDCS-R-78-931 

Performing  Organization  Name  and  Address 

Department   of   Computer   Science 

University  of   Illinois   at  Urbana-Champaign 

Urbana,    Illinois   61801 

10.  Project/Task/Work  Unit  No. 

COO-2383-0048 

11.  Contract/Grant  No. 

ENERGY/EY-76-S-02-2383 

,  Sponsoring  Organization  Name  and  Address 

Department   of   Energy 
9800   South   Cass  Avenue 
Argonne,    Illinois   60439 

13.   Type  of  Report  Si  Period 
Covered 

Research 

14. 

.  Supplementary  Notes 

.  Abstracts 

It    is   proved    that,    for   any  multistep    formul 

constant   stepsize,    there  exist   constants  b  • 

of  adjacent   steps    satisfies   b   <   h   /h      ,    <   B 

—     n     n-1  — 

of   the   formula   is   stable.      A  practical   step- 

a  which   is   strongly   stable   at 
<   1   <   B   such   that    if    the   ratio 
,    the  variable-step   implementation 
-order   control   restriction  which 

guarantees  stability  is  given. 


Key  Words  and  Document  Analysis.     17a.   Descriptors 

variable-step 

ordinary  differential   equations 

stability 

multistep  methods 

variable-order 

•  Identif iers/Open-Ended  Terms 


•  COSATI   Fie  Id  /Group 
tAvT 


Availability  Statement 

unlimited 


19.  Security  Class  (This 
Report) 

UNCLASSIFIED 


20.  Security  Class  (This 
Page 

UNCLASSIFIED 


21.   No.  of  Pages 
13 


22.   Price 


CM  NTIS-15   (  10-70) 


USCOMM-DC    40329-P7  1 


SEP  2  1  t978 


\m 


